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Abstract 

A new discrete- velocity model is presented to solve the three-dimensional Euler equations. The 
velocities in the model are of an adaptive nature — both the origin of the discrete-velocity space 
and the magnitudes of the discrete-velocities are dependent on the local flow — and are used in 
a finite volume context. The numerical implementation of the model follows the near-equilibrium 
flow method of Nadiga and Pullin [1] and results in a scheme which is second order in space (in the 
smooth regions and between first and second order at discontinuities) and second order in time. 
(The three-dimensional code is included.) For one choice of the scaling between the magnitude 
of the discrete-velocities and the local internal energy of the flow, the method reduces to a flux- 
splitting scheme based on characteristics. As a preliminary exercise, the result of the Sod shock-tube 
simulation is compared to the exact solution. 

1. Introduction 

A discrete velocity gas is an ensemble of particles with each particle taking on one of a small 
finite set of allowable velocities [2,3]. Further, the interaction between particles is defined to achieve 
the desired macro-behavior of the system, which is usually a set of partial differential equations. 
Such a discretization of the velocity space and definition of the particle interactions (collisions 
or relaxation or more generally redistribution) also form the basis for the lattice gas and lattice 
Boltzmann techniques which have been developed over the last eight years [4,5,6 and references 
therein]. In spite of the inherently compressible nature of discrete- velocity gases, lattice gases and 
lattice-Boltzmann techniques, their applications in fluid flow modeling have been restricted to the 
incompressible or very low Mach number regimes [4,5,6 and references therein, 7, 8, 9]. In this article, 
we present a discrete-velocity gas which for the first time solves the compressible Euler equations 
without introducing any artifacts of the velocity discretization over a wide range of Mach numbers. 
This is achieved by letting the discrete- velocities of the model adapt to the local flow conditions 
[10]. We show also how this scheme can be reduced to a characteristic-based flux-splitting scheme 
for the Euler equations [11,12]. 

In the next section, the model is introduced and shown to reproduce the Euler equations. 
Section 3 gives a step by step numerical evolution of the model on the lines of Nadiga and Pullin 
[1], and ?? presents a the Sod shock tube simulation [13] with the model. In ?? we show the 
reduction of the method to characteristic-based flux-splitting and end with some remarks in ??. 
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2. The Discrete- Velocity Model 

The discrete- velocity model we consider has 27 velocities (q a , a = 0, 1, ... , 26): 

q a = {q ax , q ay , q az ), q ay , q az ) e (-g, 0, q). (l) 

At any given point in space, the model thus has four different speeds 0,q,\/2q, and, \/3q. There 
is one velocity with zero speed, six velocities with speed-^, twelve velocities with speed-\/2<?, and 
eight velocities with speed-\/3<?. Note that this model consists of the more familiar two-dimensional 
nine- velocity model [14,15,16] at three different values of velocity in the z-direction ( — q,0,q), thus 
comprising 27 velocities. Collisions between particles in the model are such that they individually 
conserve mass, momentum, and energy. In particular, there are two types of collisions which 
not only conserve mass, momentum and energy but which also change the speeds of the particles 
involved. (The post-collision pair of speeds is not a mere permutation of the pre-collision speeds.) 
Describing them in the plane of the collision, the first type involves a speed- \plq particle colliding 
with a stationary particle to result in two mutually perpendicular speed-^ particles. The second 
type involves a speed-\/3<? particle colliding with a stationary particle to result again in a pair of 
particles moving mutually perpendicularly, but now one with speed q and the other with speed 
V2q. 



2.1 The Stationary Equilibrium Distribution 

Thermodynamic equilibrium of the model is defined as the state of detailed balance of all 
possible collisions. Consider the stationary equilibrium of such a gas: since there are no preferred 
directions, the particles are identified by their speeds and there are thus four variables 7io,7ii,7i2, 
and 713, where too is the probability that a particle is stationary, n\ is a sixth of the probability that 
a particle has a speed q, since there are six speed q velocities in the model, 712 is a twelth of the 
probability that a particle has speed \/2q, and 713 is an eigth of the probability that a particle has 
speed \/3<?- The subscripts in too, . . .,713 are the square of the speed divided by the unit of speed 
q. Detailed balancing of collisions results in the equilibrium condition 

n n 2 =n\, 
n n 3 =72472.2 • 

The population densities (710,711,712, and 713) are further constrained to satisfy the specified hydro- 
dynamic quantities mass 71 and energy 7ie: 



71 = no + 6toi + 12to2 + 8713, 
7ie = q 2 (3ni + 12to2 + I271.3). 



(3) 
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The stationary (n a q a = 0) equilibrium velocity distribution is obtained by the solution of above 
four equations (2) and (3): 



n f n „ e \ / e 

. 3 

n e 



713 = V \j 



(4) 



2.2 The Euler Equations 

The equations we seek to model are the Euler equations which describe the compressible inviscid 
flow of a fluid (here an ideal monatomic gas). Written in a conservative form [17], and with no 
other body forces, they are 



dp 
dt 



dpu 

~~dT~ 
dpe t 



7 ■ pu = 0, 

V • (pi + pu (8) u) = 0, (5) 



V-{(p + pe t )u} = 0. 



dt 

where ® represents the binary outer product operator, I is the unit tensor, and e t = e + u 2 /2. The 
pressure p is given by the ideal gas equation of state for a monatomic gas: p = 2/3pe. 



2.3 The Locally Adaptive Discrete Velocities 

To represent the above equations as a discrete- velocity gas, we consider the stationary (n a q a = 
0) 27-velocity gas discussed above under two simple transformations (see Fig. 1): 

• The origin of the discrete-velocity space is translated to u(x,f), where u(x,f) is the tem- 
porally and spatially varying velocity field of the Euler equations. 

• The unit of discrete velocity, q is determined locally from the specific internal energy field 
e(x,f) of the Euler equations: 

q(y.,t) = y/ae(y.,t) (6). 

The scaling factor a is a parameter in the model. To insure positivity of the distribution 
(4), the restriction on a is 

l<a<oc. (7) 
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Fig. 1 The origin of the discrete- velocity space is determined by the local flow velocity u=(u, v, w) 
of the Euler system and the unit of discrete-velocity q by its specific internal energy. The 
original nine velocities (q ax , lay, 0) where (q ax , Qay)E{ — Q, 0, q)) are shown here after they have 
adapted themselves to the local macroscopic state. The schematic has been drawn on the 
q az = w plane, where w is the z-component of the macroscopic flow velocity used in the 
Euler equations. Note that the allowable discrete velocities are now different at each point 
in space and are different at the same point with time. 

Figure 1 shows the original nine velocities (q ax , lay, ®',(<la.x, Qay)^( — 1, 0, q)) after they have adapted 
themselves to the local macroscopic state. The schematic is a projection on to the q az = w plane, 
where w is the macroscopic velocity in z-direction used in the Euler equations (5). With this kind 
of adaptation, the allowable velocities in the model are now c a (x,i) = q a (x,f) + u(x,f), i.e., the 
allowable discrete- velocities are completely different at each point in space, and are different at the 
same point in space at different times. 
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2.4 The Equivalence of the model to the Euler Equations 

The density and energy of the discrete- velocity gas are set equal to that in (5) and again note 
that n a q a has been assumed 0. The equivalence of the discrete- velocity gas system to the system 
of Euler equations might now be obvious. However, for the sake of completeness, we explain 
the equivalence. Considering the evolution of the discrete-velocity gas, the (model) Boltzmann 
equations [3,15] — a statement of the conservation of the number of particles with a particular 
dis cr et e- velo city — are 

On 

+ c a ■ Vn a = <5a(n, n), a=0,...,26, (8) 

where Q a is the nonlinear collision operator and the left hand side represents streaming of particles 

with velocity c a . The zeroth, first, and second order velocity moments of (8) give respectively, 

noting that the moments of the collision terms on the right-hand side vanish owing to the mass, 

momentum, and energy conserving nature of each collision, 

dn 

V • n a c a = 0, 



dt 



V • n a c a ® c a = 0, (9) 



V • n a -£-c a = 0, 



dt "2 

where the overbar denotes averaging with respect to the discrete- velocities. From the translation 
of the origin of the discrete- velocity gas, 



n a c a = n a (u + q a ) = pu, since n a q a = 0. 

(10) 



n a c a ® c a =n a (u + q a ) ® (u + q a ) 



-nu ® u + n a q a ® q a = pi + pu 



n a —c a =n a { h u • q a )(u + q a ) 



q 2 , u 2 



=TC a -^ u + P~^ u + n a (q a ■ u)q a 



=pe t u + n a (q a ® q a )u = {p + pe t )u. 
The equivalence is thus complete. Note that in the above equation, since n a q a = 0, the particles 



are distributed symmetrically with respect to q a and so n a q a ® q a is isotropic and reduces to the 
pressure p. 

For convenience, the Euler equations (5) and the moment equations (9) may be rewritten as 

df 




with f = ( pu ] , (11) 
G = [pu, pi + pu ® u, {p + pe t )u] 



3. The Numerical Technique 



Hyperbolic systems of conservation laws like the Euler equations (5) in the present case, admit 
weak solutions in the form of shocks. At these shocks, the gradients of the primary quantities 
are infinite, and obviously this cannot be represented correctly in a shock-capturing numerical 
scheme. Shock-capturing schemes are ones in which all grid points are treated in exactly the 
same fashion irrespective of whether they are inside a shock or outside, and as opposed to shock- 
tracking schemes which keep track of where the shocks are and treat grid points located in shocks 
differently from the others. While each have their advantages and disadvantages, used on present 
day supercomputers, shock-capturing methods are clearly preferable because of their homogenity of 
computation. Numerically capturing shocks in nondissipative systems like the Euler equations poses 
the problem of dispersive ripples: a phenomenon in which as a wave steepens, energy flows into the 
smaller scales and since there is no dissipation, accumulates in the smallest allowed wavelengths — 
those of the grid-spacing. Thus with the formation of any shock, energy piles up in grid scale 
oscillations and swamps out scales of interest. A way out of this problem is to dissipate energy 
at the smallest length scales — the grid-scales — either explicitly or implicitly. Explicit artificial 
viscosity in general needs adjustments depending on the problem, while physics-based (as opposed 
to the inviscid idealization represented by the Euler equations) implicit artificial (from the point of 
view of the Euler equations) viscosity is robust and results in numerical shock widths of the same 
order as the actual viscous (Navier-Stokes in the present case) shock widths. 

We use exactly such an implicit artificial viscosity technique (also called a kinetic numerical 
scheme owing to the physical kinetic basis of the scheme) which was developed in Nadiga and Pullin 
[1] for the discrete- velocity framework. A brief description of its usage here follows: For simplicity, 
consider the computational domain divided into uniform cubical cells, at the centroids of which 
are stored the cell-averaged values of (p,pu, and pe t ). The evolution at each centroid proceeds as 
follows: 

Step 1: Calculate the local unit of discrete- velocity using (6) and (7). 

Step 2: Using (4) and the macroscopic variables of density p and specific internal energy e, 
calculate the population densities of the four speeds no,ni,n2, and 72.3 . Note that in light of 
restriction (7) and owing to the fact that the stationary distribution function was calculated based 
on detailed balancing, the above four population densities are necessarily positive. 

Step 3: Calculate the split fluxes G + and G~ at the centroids using the definitions 



- 8- 



where n°a c ax is the average n a c ax taken only over the discrete velocities c a which have a positive 
ce- component c ax , etc For example, when u > 0, and w > 0, but v < 0, 



G 



^uP + (u + q)Q u 2 P + (u + q) 2 Q uvP + (u + q)vQ uwP + (u + q)wQ n c a ax>0 ^c ax ^ 

(v + q)Q (v + q)uQ (v + q) 2 Q (v + q)wQ n a ay> °fc ay > 

\wP + (w + q)Q wuP + (w + q)uQ wvP + (w + q)vQ w 2 P+(w + q) 2 Q n c a az>0 ^fc az / 

(13) 



G 



(u-q)Q {u-q) 2 Q {u-q)vQ (u - q)wQ n c a ax<0 fc ax ^ 

vP+(v-q)Q vuP + (v-q)uQ v 2 P+(v-q) 2 Q vwP + (v - q)wQ n a ay< °^fc ay 

\ {w-q)Q {w-q)uQ {w-q)vQ {w - q) 2 Q n c a az<0 ^fc az / 

where P = (no + 4ni + 4^) is the density of particles in say the c ax = u plane and Q = {n\ + 
4^2 + 4713) is the density of particles in say the c ax = u + q or c ax = u — q plane. In the above two 

2 2 

equations, ria ax>0 ^fc ax , . . ., ria az<0 ^fc az have been left as such for compactness of notation. 

Step 4: Assuming a linear distribution of the fluxes within the cells in the direction under consid- 
eration, interpolate the split fluxes G + and G~ to the cell boundaries: 

• interpolate Gf x , G^ 2 , G^ 3 , G^ 4 , and G^ 5 to (i + \,j,k), 

• interpolate G 21 , G 22 , G 23 , G 24 , and G 25 to — \,k), and so on. (Note that the first sub- 
script of G corresponds to the coordinate direction and the second to one of the conserved 
quantity.) 

and apply the minmod limiter to the interpolated fluxes: 

G n( i + \ii' k ) = G ti( i >3> k ) + ^mmmod(A bck G^ 1 (i,j,k),A fwd P^ 1 (i,j,k)), 

G 2i(hJ ~ \i k ) = G 21 (i,j,k)- ^mmmod(A bck G 21 (i,j,k),A fwd G 21 (i,j,k)), (14) 
and so on, and where 

AfvdG+^iJ^) = G+ 1 {i+l,j,k)-G+ 1 {i,j,k) 

is the first forward difference of in the cc-direction at the centroid of cell k). 

A bck G 21 (i,j,k) = G 21 (i,j,k) - G 21 (i,j - 1, k) 

is the first backward difference of G 21 in the y-direction at the centroid of cell (i,j,k) and so on. 
And minmod is the one dimensional total- variation-diminishing operator as discussed in [18,19]: 
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with sgn(p) being the sign of p and \p\ the absolute value of p. 
Step 5: Calculate the fluxes at the cell boundaries: 

G(i + i) = G+(i + i) + G-(i + i) 

where (i + \) = {i + k), + \,k), and k + \) in turn. 

Step 6: Advance the primary variables by one half of the time step to get the midpoint values: 

fto+At/2 = ft + ^ ( V . G*°) , (16) 

, „ (^\ G 11 (i + lj,k)-G 11 (i-lj,k) t 
where e.g. V • G21 = 7 h 

G 2 i{h3 + h k )~ G 2i{h3 ~ \,k) G 3 i(i,j,k+ \) - G 3 i(i,j,k- \) 
Ay Az 

Step 7: Repeat steps 1-4 using the midpoint values f'o+At/2 ^ Q obtain the fluxes at the midpoint 
G'° +At / 2 , and take the full time step to obt ain the new time level values f' 1 : 

f*i = fh + Ath- G t0+At/2 ) . (17) 



In the above procedure the fluxes were interpolated and limited, instead, the primary quantities 
f could be interpolated and limited. We have done both and the behavior of the two are essentially 
identical. The time step in the method is limited by the CFL stability criterion 

max ± l\-^j^J < !> U £ (u,v,w) and AX £ (Ace, Ay, Az). (18) 

The code (in Connection Machine Fortran) described by the above step-by-step procedure is in- 
cluded in the Appendix. Notwithstanding the somewhat complicated procedural description above, 
the simplicity of the resulting code is evident. 



4. A Numerical Example 

The code presented in the Appendix, a Connection Machine Fortran implementation of the 
second order (in space and time) scheme described in the previous section, was used to calculate 
the Sod shock-tube problem. The Sod test case has the initial conditions (pi = 1, ui = 0,pi = 1) and 
(p r = 0.125,1^ = 0,p r = 0.1) corresponding to an initial pressure ratio of 10 and a density ratio of 
8. Subscript I denotes the left half and subscript r denotes the right half at time 0. For a monatomic 
gas, this reduces to the initial conditions (pi = 1, ui = 0, e/ = 1.5) and (p r = 0.125, u r = 0, e r = 1.2). 
Only 128 points were used for this simulation, at a timestep corresponding to a CFL number of 0.69. 
(The a in (6) was set at 10/9.) In Fig. 2, the exact density, specific internal energy, velocity and 
pressure profiles (solid lines) are compared to the corresponding profiles (open diamonds) obtained 
from the discrete velocity simulation. The agreement is good: The shock is typically 3-4 cell-widths 
and the edges of the rarefaction are not too badly rounded. The spreading of the contact surface 
is however substantial as with other flux-splitting schemes. 
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Fig.2 The Sod shock-tube problem: At a dimensionless time of 0.288 after the diaphragm burst 
(the diaphragm is initially at x=0), the inital discontinuity in pressure and density (pressure 
ratio of 10 and a density ratio of 8) is resolved into a right-going shock wave located at about 
x=0.5 , a left-going rarefaction centered about x=-0.3, and a right-going contact surface at 
about x=0.2. The solid line is the exact solution and the diamonds are the result of the 
simulation using the locally adaptive discrete velocities. 

5. The recovery of Characteristic-based fiux-splitting 

From the point of view of finite differencing, high-order shock-capturing techniques used for 
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numerically solving inviscid hyperbolic systems of conservation laws are in general based on high- 
order upwinding. Physically, upwinding can be viewed in two different ways: Firstly, as resulting 
from a Gudonov methodology of resolving discontinuities in data at cell interfaces using an approx- 
imate Riemann solver [11,12]. And secondly, as resulting from the local solution of the collisionless 
Boltzmann equations [20,21,22]. The present scheme falls into the second category but differs from 
other such schemes in having a discretized velocity space. 

If the unit of discrete- velocity <?(x,i) is set equal to the local sound speed a(x,i), 



where 7 = 5/3 for a monatomic gas. Next, from the formulae for the split fluxes (13) with P = OAp, 
Q = 0.3/O, and q = a, in the one dimensional case, there are three beams at the three characteristic 
speeds u + a, u, and u — a. Now, depending on the signs of u + a, u, and u — a, the positive fluxes 
consist of the beams whose speeds are positive and the negative fluxes consist of the beams whose 
speeds are negative. Thus, the scheme now exactly propagates information along the characteristics 
of the Euler equations and reduces to a characteristic-based flux-splitting scheme [11,12]. In higher 
dimensions however, the inherent multi- dimensional (and upwinding) nature of particle motion in 
the present scheme appears to make it different and needs to be further investigated. 



6. Conclusion 

Wanting to incorporate the simple and elegant way in which discrete-velocity gases (which 
include the lattice gas and lattice-Boltzmann formulations) combine physics and numerics in the 
finite-volume techniques for the Euler equations led us to an adaptive discrete velocity model for 
the Euler equations. The adaptive nature of the discrete velocities (i.e. the variable origin and 
local scaling of the discrete velocities) seems to bridge the gap between the newer discrete- velocity 
techniques and the more conventional flux-splitting techniques. Based on earlier work, we develop 
a second order (in space and time) scheme which is very simple and yet robust: it can handle 
flows over a wide range of Mach numbers accurately and capture shock jumps over 3-4 slab widths 
with no oscillations. We are presently working on an extension of the method wherein there is a 
Bhatnagar-Gross-Krook [23] like relaxation of the particle distributions over the timestep At, on 
the lines of [22]. This will further reduce the diffusive nature of the shocks and the rounding of the 
corners of rarefaction waves. We wish to emphasis at this point that the new scheme being a finite 
volume technique, the use of discrete-velocity gases is now possible with arbitrary and irregular 
spatial meshes. (Though this point was implicit in [1,10], it has sometimes been overlooked.) The 
limiting scheme used to achieve second order accuracy is however one- dimensional and therefore the 
second order scheme suffers from all the drawbacks of dimension- split methods. It remains to be 
investigated if the method can be made genuinely multi-dimensional. We also plan to investigate 
the extension of the model to solve the Navier-Stokes equations from its present capability of solving 
the Euler equations. 




(19) 
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The full three-dimensional code is presented in the Appendix. It is in CMFortran (a CM 
extension of Fortran90) for ease of understanding and presentation. A Fortran77 version of the 
code may be had from the author upon request. This work was supported by DOE in part under 
the CHAMMP program. 
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THE CMF CODE LISTING 



THE MAIN ROUTINE 

include ' /usr/ include/ cm/ cmssl-cmf .h' 

include ' include . h ' 

real ,array(5 ,nx,ny,nz) : :fnew,fold 
cmf$ layout fnew( :serial, , ,) ,fold( :serial, , ,) 
c 

c INITIALIZE HERE. (Code left out.) 
c 

call CMF_cm_array_to_f ile(22 ,f new , ioerr) 
cf lav=0 
nav=0 
print * , ' ' 
c MAIN LOOP TO DO MIDPOINT INTEGRATION 
do 100 kt=l,nt 

call flux3d(fnew) 

f old=f new+0 . 5*dt*df dt 

call flux3d(fold) 

f new=f new+dt*df dt 

if (mod(kt ,ntwrite) .eq.0) then 
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print 300,kt*dt ,kt , cf lav/nav,cf 1 
call CMF_cm_array_to_f ile(22 ,f new, ioerr) 
endif 
100 continue 

300 format ('time=' ,f8.4, ' kt=',i5,' cf lav= ' ,f 7 . 3 , ' cfl=',f7.3) 
stop 
end 



SUBROUTINE FLUX3D 

c This subroutine calls fluxld 3 times, once for each dimension. 

subroutine flux3d(f basic) 

include " include. h" 

real ,array(5 ,nx,ny,nz) : :fbasic 
cmf$ layout f basic ( : serial, , , ) 

dfdt=0.0 

call f luxld(2,3,4,fbasic) 
call f luxld(3,4,2,fbasic) 
call f luxld(4,2,3,fbasic) 
return 
end 



SUBROUTINE FLUX1D 

c For the cell centered at i, this subroutine calculates the 

c forward going fluxes at i+1/2 and backward going fluxes at i-1/2 

c The limiting procedure is used on the primary quantities. 

subroutine f luxld(iaxis , j ax,kax,f basic) 

include " include. h" 

integer iaxis, jax,kax 

real,array(nx,ny,nz) : :tmpll ,tmpl2,tmpl3 
real,array(nx,ny,nz) : :rho,u,v,w,e,q 
real,array(5,nx,ny,nz) : :fbasic,gxp,gxm,dff 
cmf$ layout f basic ( : serial, , ,) 

cmf$ layout gxp( : serial ,,,) ,gxm( : serial ,,,) ,dff (: serial ,,, ) 

cmf$ layout tmpll ( , , ) ,tmpl2( , , ) ,tmpl3( , , ) 

cmf$ layout rho( , , ) ,u( , , ) , v( , , ) ,w( , , ) , e( , , ) ,q( , , ) 
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call minmod(f basic , iaxis ,dff ) 
gxp=f basic+O . 5*df f 
gxm=f basic-0 . 5*df f 

Calculate Forward going fluxes at i+1/2 (ru ,ruu ,rvu ,rwu ,retu) 



rho=gxp ( 1 , : 
u=gxp( iaxis 
v=gxp(jax, 
w=gxp(kax , 
e=gxp(5, : , 



) 



: ) /rho 
: ) /rho 
: ) /rho 
: ) /rho 
e=e-0 . 5* (u*u+v*v+w*w) 
q=sqrt (alpha*e) 
c Note: correct only for ID 

cf l=maxval (max (abs (u-q) , abs (u+q) ) ) *dt/dxyz (iaxis) 
cf lav=cf lav+cf 1 
nav=nav+l 
e=e/(q**2+l .Oe-30) 
c Calculate the stationary equilibrium distribution at i+1/2 
tmp=3 . -2 . *e 
rho=rho/27 . 
n0=rho*tmp**3 . 
nl=rho*tmp*tmp*e 
n2=rho*tmp*e*e 
n3=rho*e**3 . 



:)=4.*(n2+n3)+nl 
:)=4.*(n2+nl)+n0 



c Mass 

tmp 1 ( 1 , : : 
tmp2 ( 1 , : : 
c X-Momentum 

tmp 1 (iaxis )=tmpl (1 , :,:,:)* (u+q) 
tmp2( iaxis )=tmp2(l , :,:,:)* (u) 
c Y-Momentum 

tmpl ( j ax , : , : , : )=tmpl (1 , :,:,:) *v 
tmp2( j ax , : , : , : )=tmp2(l , :,:,:) *v 
c Z-Momentum 

tmpl (kax , : , : , : )=tmpl (1 , :,:,:) *w 
tmp2 (kax ,:,:,:) =tmp2 ( 1 ,:,:,:) *w 

c Energy 

tmpll=v*v 
tmpl2=w*w 
tmpl3=2.*q*q 

tmp=2.*(2.*n3+n2)*(tmpll+tmpl2+tmpl3)+(2.*n2+nl)*(tmpll+tmpl2) 
tmpl (5, : , : , : )=0 . 5* (tmpl (1 , : , : , : ) *( (u+q) *(u+q)+tmpll+tmpl2)+tmp) 
tmp2(5, : , : , : )=0 . 5* (tmp2(l , : , : , : )*(u*u+tmpll+tmpl2)+ 
> 2.*(2.*n2+nl)*(tmpll+tmpl2+tmpl3)+(2.*nl+n0)*(tmpll+tmpl2)) 



tmpl2=(u+abs(u)) 
forall(i=l:5) 

gxp(i, : , : , :)=(u+q)*tmpl(i, : , : , :)+0.5*(tmpl2*tmp2(i, :,:,:)) 
Calculate Backward going fluxes at i-1/2 (ru ,ruu ,rvu ,rwu ,retu) 
rho=gxm(l, :,:,:) 
u=gxm(iaxis ) /rho 
v=gxm(jax, : , : , : )/rho 
w=gxm(kax, : , : , : )/rho 
e=gxm(5 ) /rho 
e=e-0 . 5* (u*u+v*v+w*w) 
q=sqrt (alpha*e) 
e=e/(q**2+l .Oe-30) 

cf l=maxval (max (abs (u-q) , abs (u+q) ) ) *dt/dxyz ( iaxis) 

cf lav=cf lav+cf 1 

nav=nav+l 

tmp=3 . -2 . *e 

rho=rho/27 . 

n0=rho*tmp**3 . 

nl=rho*tmp*tmp*e 

n2 =rho * tmp * e * e 

n3=rho*e**3 . 



tmpl(l , : , : , 


) = 


=4.*(n2+n3)+nl 




tmp2 (1,:,:, 


) = 


=4.*(n2+nl)+n0 




tmp 1 (iaxis , 




, :)=tmpl(l, 




, :)*(u-q) 


tmp2( iaxis , 




, : ) =tmp2 ( 1 , 




,:)*(u) 


tmpKjax, : , 




)=tmpl(l, : , 




)*v 


tmp2(jax, : , 




) =tmp2 ( 1 , : , 




)*v 


tmpl(kax, : , 




)=tmpl(l, : , 




) *w 


tmp2(kax, : , 




) =tmp2 ( 1 , : , 




) *w 


tmpll=v*v 








tmpl2=w*w 








tmpl3=2.*q*q 







tmp=2.*(2.*n3+n2)*(tmpll+tmpl2+tmpl3)+(2.*n2+nl)*(tmpll+tmpl2) 
tmpl(5, : , : , : )=0 . 5* (tmpl (1 , : , : , : ) *( (u-q) *(u-q)+tmpll+tmpl2)+tmp) 
tmp2(5, : , : , : )=0 . 5* (tmp2(l , : , : , : )*(u*u+tmpll+tmpl2)+ 

2.*(2.*n2+nl)*(tmpll+tmpl2+tmpl3)+(2.*nl+n0)*(tmpll+tmpl2)) 
tmpl2=(u-abs (u) ) 
forall(i=l:5) 

gxm(i , : , : , : )=(u-q) *tmpl (i , : , : , : )+0 . 5* (tmpl2*tmp2(i , :,:,:)) 
tmpl=cshif t (gxm, iaxis , 1) 
tmp2=cshif t (gxp , iaxis , -1) 
df dt=df dt+ (tmp2-gxp+gxm-tmpl) / dxyz (iaxis) 
return 



- 17- 



end 



SUBROUTINE MINMOD 

c Given an array y, this routine returns the limited first difference 
c using the minmod limiter. 

subroutine minmod (y ,j ,ydff) 

include " include. h" 

real ,array(5 ,nx,ny,nz) : :y,ydff 
cmf$ layout y (: serial ,,,) ,ydff (: serial ) 

data eps/1.0e-30/ 

tmpl=(cshift(y, j , l)-y) 

tmp2=(y-cshift(y, j ,-1)) 

ydf f =tmp2*tmpl 

ydf f =sign(l . , turpi) *max(0 . ,ydf f ) / abs (ydf f +eps) 
ydff=ydff *min( (abs (turpi) ) , (abs(tmp2))) 
return 
end 



INCLUDE. H 

c Note that the number of points in the x direction, nx is 

c set to 256 since we are running with periodic boundary conds . 

c Only half of that i.e. 128 points are really used. 

integer .parameter : :nx=256 ,ny=16 ,nz=l 

real , array (nx,ny,nz) : :nO,nl J n2,n3,tmp 

real ,array(5 ,nx,ny,nz) : : df dt ,tmpl ,tmp2 
cmf$ layout n0( , , ) ,nl ( , , ) ,n2( , , ) ,n3( , , ) ,tmp( , , ) 
cmf$ layout dfdt( :serial, , ,) ,tmpl( :serial, , ,) ,tmp2( :serial, , ,) 

common/ scalar/ dt , qurat , alpha , cf 1 , cf lav ,nav , dxyz (4) 

common/vector/ nO ,nl ,n2 ,n3 ,tmp ,df dt , turpi ,tmp2 



